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Abstract 

Continuous-time  birth-death  Markov  proeesses  serve  as  useful  models  in  population  biol¬ 
ogy.  When  the  birth-death  rates  are  nonlinear,  the  time  evolution  of  the  first  n  order  mo¬ 
ments  of  the  population  is  not  elosed,  in  the  sense  that  it  depends  on  moments  of  order 
higher  than  n.  For  analysis  purpose,  the  time  evolution  of  the  first  n  order  moments  is  often 
made  to  be  elosed  by  approximating  these  higher  order  moments  as  a  nonlinear  funetion  of 
moments  up  to  order  n,  whieh  we  refer  to  as  the  moment  closure  function. 

In  this  paper,  a  systematie  proeedure  for  eonstrueting  moment  elosure  funetions  of  ar¬ 
bitrary  order  is  presented  for  the  stoehastie  logistie  model.  We  obtain  the  moment  elosure 
funetion  by  first  assuming  a  eertain  separable  form  for  it,  and  then  matehing  time  deriva¬ 
tives  of  the  exaet  (not  elosed)  moment  equations  with  that  of  the  approximate  (elosed) 
equations  for  some  initial  time  and  set  of  initial  eonditions.  The  separable  strueture  ensures 
that  the  steady-state  solutions  for  the  approximate  equations  are  unique,  positive  and  real, 
while  the  derivative  matehing  guarantees  a  good  approximation,  at-least  loeally  in  time. 
Moreover,  the  aeeuraey  of  the  approximation  ean  be  improved  by  inereasing  the  order  of 
the  approximate  model.To  the  best  of  our  knowledge,  this  paper  is  the  first  to  propose  a  sys¬ 
tematie  proeedure  to  eonstruet  moment  elosure  funetions  of  arbitrary  order  that  guarantee 
biologieally  meaningful  equilibria. 

A  host  of  other  moment  elosure  funetions  previously  proposed  in  the  literature  are  also 
investigated.  Among  these  we  show  that  only  the  ones  that  aehieve  derivative  matehing 
provide  a  elose  approximation  to  the  exaet  solution.  Moreover,  we  improve  the  aeeuraey 
of  several  previously  proposed  moment  elosure  funetions  by  foreing  derivative  matehing. 
However,  for  eertain  ranges  of  parameter  models,  moment  elosure  funetions  that  laek  the 
separability  property  lead  to  biologieally  meaningless  seenarios  of  imaginary  and  even  sta¬ 
ble  negative  steady-states  of  the  elosed  moment  equations. 
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1  INTRODUCTION 


Continuous-time  birth-death  Markov  processes  have  been  extensively  used  for  mod¬ 
eling  stochasticity  in  population  biology  (Matis  and  Kiffe,  2002,  1996;  Matis  et  ah, 
1998).  The  time  evolution  of  these  processes  is  typically  described  by  a  single  equa¬ 
tion  for  a  grand  probability  function,  where  time  and  species  populations  appear  as 
independent  variables,  called  the  Master  or  Kolmogorov  equation  (Bailey,  1964). 
However,  this  equation  can  only  be  solved  for  relatively  few,  highly  idealized  cases 
and  a  more  reasonable  goal  is  to  determine  the  time  evolution  of  a  few  low-order 
statistical  moments. 


In  this  paper,  a  method  for  estimating  low-order  statistical  moments  is  introduced 
for  the  stochastic  logistic  model.  This  model  was  first  introduced  by  Barlett  et  al. 
(1960)  and  is  a  continuous-time  birth-death  Markov  process  involving  a  single 
specie,  with  birth  and  death  rates  being  polynomials  of  order  2.  Although,  one  can 
directly  use  the  Kolmogorov  equation  to  derive  differential  equations  for  the  time 
evolution  of  moments  of  the  process,  in  this  paper  we  use  an  alternative  method. 
We  model  the  stochastic  logistic  model  as  a  Stochastic  Hybrid  System  (SHS)  whose 
state  X  is  the  population  of  the  specie.  Then,  the  time  evolution  of  the  moments  is 
obtained  using  results  from  the  SHS  literature  (e.g.,  Hespanha,  2004).  Details  of 
the  stochastic  logistic  model  and  its  modeling  as  a  SHS  are  presented  in  Section  2. 


Table  1 

Separable  derivative-matehing  moment  elosure  funetion  <p*+i(Ai)  for  n  G  {2,3,4}. 


n=2  n=3  n=4 


3“ 


<Pn+l(Ai)  § 

_ 


Let  Hin  be  the  m'*  order  uncentered  moment  of  x,  i.e.  limit)  =  E[x(t)'"]  where  x(t) 
denotes  the  population  of  the  specie  at  time  t.  We  will  show  in  Section  3  that  for 
the  stochastic  logistic  model  the  time  derivative  of  flm  is  a  linear  combination  of 
moments  up  to  order  m  -f  1 .  Hence,  if  one  stacks  all  moments  in  an  infinite  vector 
/ioo  =  [ill , /i2,  ■  ■  ■  ]^>  its  dynamics  can  be  written  as 

floo  -  AcxDjttoO,  (1) 

for  some  infinite  matrix  Aoc.  As  the  above  infinite  dimensional  system  cannot  be 
solved  analytically,  we  truncate  (1)  by  creating  a  vector  /i  =  [/ii, . . .  where  n 

is  the  order  of  the  truncation.  Its  dynamics,  given  by 


li=Aii  +  Bii„+i  (2) 

for  some  matrices  A  and  B  is  not  closed,  in  the  sense  that  the  time  evolution  of 
the  vector  /i  depends  on  the  n  +  order  moment  /i„+i.  For  analysis  purposes,  we 
close  the  above  system  by  approximating  /i„+i  as  a  nonlinear  function  (p„_|_i(/i)  of 
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moments  up  to  order  n.  This  procedure  is  commonly  referred  to  as  moment  closure. 
We  call  as  the  moment  closure  function  for  Let  v  = 

denote  the  state  of  the  new  closed  system.  Then,  its  dynamics  is  given  by 

v=Av  +  5(p„+i(v)  (3) 

and  is  referred  to  as  the  truncated  moment  dynamics.  We  denote  the  states  of  (2) 
and  (3)  using  different  symbols  because  ji  refers  to  the  actual  moment  dynamics 
whereas  V  to  an  approximated  moment  dynamics. 

In  Section  4,  we  consider  moment  closure  functions  which  have  the  following  sep¬ 
arable  form: 


for  appropriately  chosen  constants  G  M.  These  constants  are  obtained  by  match¬ 
ing  time  derivatives  of  /i„+i  and  (/i)  in  (2)  and  (3)  respectively,  at  some  initial 
time  to,  for  a  deterministic  initial  condition  x(to)  =  xq  with  probability  one.  The 
reason  for  this  lies  in  the  fact  that  the  class  of  deterministic  distributions  forms  a 
natural  basis  for  the  infinite  dimensional  space  containing  the  vector  /ioo.  We  refer 
to  this  moment  closure  as  the  separable  derivative -matching  moment  closure.  We 
show  that  for  all  n  >2,  this  determines  the  function  uniquely,  which  is  in¬ 
dependent  of  the  birth  and  death  rates.  Table  1  shows  the  functions  that  we 
obtained  for  truncations  of  order  n  =  2,3  and  4.  The  striking  feature  of  the  sepa¬ 
rable  derivative-matching  moment  closure  is  that  the  accuracy  of  the  approximate 
moment  dynamics  improves  by  increasing  order  of  truncation  and  the  dependence 
of  higher  order  moments  on  lower  order  ones  is  consistent  with  x(t)  being  lognor- 
mally  distributed,  in  spite  of  the  fact  that  the  derivative  matching  procedure  used  to 
construct  did  not  make  any  assumption  on  the  distribution  of  the  population. 

Alternative  moment  closure  methods  which  have  appeared  in  literature  typically 
construct  the  moment  closure  functions  (p  by  directly  assuming  the  probability  dis¬ 
tribution  to  be  normal  (Whittle,  1957),  lognormal  (Keeling,  2000),  poisson  or  bino¬ 
mial  (Nasell,  2003a).  We  refer  to  them  as  normal,  lognormal,  poisson  and  binomial 
moment  closures  respectively  and  review  them  in  Section  5.  In  Section  6,  they  are 
compared  with  the  separable  derivative-matching  moment  closure.  In  Section  6.1 
the  comparisons  are  done  based  on  how  well  the  moment  closure  function  (jl) 
approximates  /i„+i.  Towards  that  end,  we  introduce  the  error 

°°  (t  —  t  Y 

en+\{t)  :=At„+i(t)-(p„+i(At(0)  =  E  ^  £«+i(^o) 

1=0  *• 

where  we  expanded  the  error  as  a  Taylor  series  with  e^^^(.ro)  defined  to  be 

„+lW)-  \t=tQ  \t=tQ 
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We  call  (a:o)  the  order  derivative  matching  error.  Ideally  one  would  like  this 
error  to  be  zero  but  this  is  generally  not  possible.  When  x(?o)  =  xq  with  probability 
one,  the  derivative  matching  error  is  typically  a  polynomial  in  xq.  The  lesser  the 
order  of  this  polynomial,  the  lesser  is  the  error  en+i{t),  and  hence  the  better  is 
in  approximating 


We  show  that  for  n  =  2,  all  the  above  moment  closure  functions  perform  derivative 
matching  except  the  poisson  moment  closure  function  by  Nasell  (2003a).  This  is 
because,  it  has  a  O'*  order  derivative  matching  error  £3(^:0)  which  grows  linearly 
with  Xq  while  for  separable  derivative-matching,  lognormal,  binomial  and  normal 
moment  closure  functions  the  O'*  order  error  is  always  zero.  Hence,  Nasell’s  pois¬ 
son  moment  closure  function  exhibits  a  larger  initial  error  than  the  others.  We  pro¬ 
pose  an  alternative  poisson  moment  closure  function,  for  which  £3  (a:o)  =  0,  and 
show  that  it  performs  better  than  the  one  proposed  by  Nasell  (2003a). 


Although  the  above  moment  closures  provide  good  estimates  for  a  second  order  of 
truncation  (n  =  2),  it  is  typically  beneficial  to  consider  n>  3  because  they  lead  to 
better  moment  approximations  and  reduce  the  errors  by  a  few  orders  of  magnitude. 
However,  for  distributions  like  the  normal  and  the  lognormal,  which  are  charac¬ 
terized  by  less  than  2  parameters,  n  >  3  typically  leads  to  multiple  normal  and 
lognormal  moment  closure  functions.  In  Section  6. 1 .2  we  illustrate  how  derivative 
matching  can  be  used  as  a  tool  for  gauging  performances  of  these  multiple  moment 
closure  functions  and  choosing  the  ones  among  them  that  yield  the  least  derivative 
matching  errors. 


In  particular,  we  show  that  among  the  moment  closure  functions  consistent  with  the 
lognormal  distribution,  the  separable  derivative-matching  moment  closure  function 
yields  the  least  order  polynomial  for  the  derivative  matching  error,  and  hence,  ex¬ 
hibits  the  best  performance.  Based  on  derivative  matching  we  also  propose  families 
of  normal  moment  closure  functions  which  provide  good  approximations  for 
Towards  the  end,  for  n  =  3,  we  propose  a  new  moment  closure  function,  for  which, 
unlike  the  lognormal  or  normal  moment  closure  functions  both  the  O'*  and  H'  order 
derivative  matching  error  are  zero,  and  hence,  provides  the  best  estimates  for  /i4  as 
compared  to  them,  at  least  locally  in  time. 


In  Section  6.2  comparisons  are  done  based  on  the  steady-state  solutions  of  the  trun¬ 
cated  moment  dynamics  (3).  We  show  that  the  separable  derivative  matching  mo¬ 
ment  closure,  always  yields  a  unique  non-trivial  positive  real  steady-state  Vn  G  N>2. 
Thus,  it  is  preferable  to  the  other  moment  closures  techniques  which  lack  the  sep¬ 
arable  structure  and  exhibit  spurious,  imaginary  and  even  stable  negative  steady- 
states,  which  would  be  biologically  meaningless. 
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2  Stochastic  Logistic  Model 


2.1  Model  formulation 


The  stochastic  logistic  model  is  the  stochastic  birth-death  analogous  model  of  the 
well-known  deterministic  Verhulst-Pearl  equations  (Pielou,  1977)  and  has  been  ex¬ 
tensively  used  for  modeling  stochasticity  in  population  biology  (Matis  and  Kiffe, 
2002,  1996;  Matis  et  ah,  1998).  For  this  continuous-time  birth-death  Markov  pro¬ 
cess,  the  conditional  probabilities  of  a  unit  increase  and  decrease,  respectively,  in 
an  “infinitesimal”  time  interval  fd  +  dt]  is  given  by 


P{x(t  -\-dt)  =  .r  -I-  1  |x(t)  =x}  = 


ri(x)dt,  \/x<U 


0,  otherwise 

P{x(t  -\-dt)  =x—\  |x(t)  =x}  =  x{x)dt, 

where  x(t)  G  N  represents  the  population  size  at  time  t, 

ri(x)  :=  aix  —  bix^  >  0,  x(x)  ■=  ii2X  +  b2X^  >  0,  V.rG(0,  C/) 


(4) 


and 


f/ :=  ai/Z?i  G  ai  >  0,  ^2  >  0,  bi>0,  b2>0. 

We  assume  that  the  initial  condition  satisfies  x(to)  G  { 1, 2, . . . ,  t/},  and  hence,  x(t)  G 
{0, 1, . . . ,  t/},  Vt  G  [0,oo)  with  probability  one.  We  call  U  the  population  limit. 


2.2  Stationary  and  Quasi-Stationary  Distributions 


Since  the  birth  and  death  rates  are  zero  for  x  =  0  {r\{0)  =  ;t(0)  =  0)  we  have  that 
X  =  0  is  an  absorbing  state  and  eventual  convergence  to  the  origin  is  certain.  Thus, 
the  stationary  distribution  is  degenerate  with  probability  one  at  the  origin.  However, 
though  out  this  paper  we  assume  the  mean  time  to  extinction  to  be  very  large, 
in  which  case  there  exists  a  “quasi-stationary”  distribution  (Nasell,  2001).  Barlett 
et  al.  (1960)  shows  that  a  good  approximate  for  Px  (the  probability  that  x  =  .r  at 
“quasi-equilibrium”)  can  be  numerically  obtained  using  the  following  recurrence 
relationship 


X{x)Px  =  ri{x-l)Px-i,  Vxe{2,3,...,U}.  (5) 
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2.3  Transient  Distributions 


Ideally  one  would  like  to  determine  the  exact  probability  distribution  of  x(t)  at  any 
time  t.  It  can  be  shown  that  the  probability  Px{t)  that  x(t)  =  x  satisfies  the  following 
differential  equations 


pQ{t)=x{m{t)  (6a) 

Pi{t) = ximit)  -  [T](i) +z(i)]A(o  (6b) 

Px{t)  =  x{x+  ^)Px+iit)  -  [ri{x)+xix)]Px{t)  +  ri{x-l)Px-iit),  (6c) 
Pu{t)  =  -X{U)Pu{t)  +  ri{U-  l)Pu-i{t)  (6d) 


(Bailey,  1964).  These  differential  equations  are  commonly  referred  to  as  the  Master 
or  Kolmogorov  equations.  If  the  population  limit  U  is  small,  the  above  system  of 
equations  can  be  solved  numerically.  However,  for  large  U ,  a  more  reasonable  goal 
(and  one  that  is  of  primary  interest  in  applications)  is  to  determine  the  evolution  of 
the  some  lower-order  moments  of  x(t). 


3  Time  Evolution  of  Moments 


3. 1  Modeling  the  Stochastic  Logistic  Model 


To  model  the  time  evolution  of  x(t),  we  consider  a  special  class  of  systems  known 
as  Stochastic  Hybrid  Systems  (SHS).  These  systems  were  introduced  by  Hespanha 
and  Singh  (2005)  to  model  the  stochastic  time  evolution  of  the  populations  of  differ¬ 
ent  species  involved  in  a  chemical  reaction.  More  specifically,  to  fit  the  framework 
of  our  problem,  these  system  are  characterized  by  two  reset  maps: 

X  I— >  ^i(x)  :=  x-|- 1,  X  I— >•  ^2(x)  :=  X  —  1  (7) 

one  corresponding  to  a  birth  and  the  other  to  a  death,  with  associated  transition 
intensities  given  by 


Ai(x)  :=r](x),  A2(x):=^(x).  (8) 

Between  births  and  deaths  the  population  remains  constant  and  thus  x  =  0.  In 
essence,  whenever  a  “birth  event”  or  a  “death  event”  takes  place,  the  corresponding 
reset  ^;(x)  is  “activated”  and  x  is  reset  accordingly,  furthermore,  the  probability  of 
the  activation  taking  place  in  an  “infinitesimal”  time  interval  (t,  t  +  dt]  is  determined 
by  the  associated  transition  intensities  Xi{x)dt. 
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3.2  Moment  Dynamics 


Given  m  G  {1,2, . . .},  we  define  the  order  (uncentered)  moment  of  x  to  be 


:=  E[x(t)n,  'it  >  0. 

X=1 


(9) 


The  time  evolution  of  moments  is  given  by  the  following  result,  which  is  a  straight¬ 
forward  application  of  Theorem  1  in  Hespanha  (2004)  to  the  above  SHS. 

Theorem  1:  The  time  evolution  of  jim  is  given  by 


dflfr, 

dt 


E 


1=1 


(10) 


Using  the  above  Theorem,  we  show  in  Appendix  A  that  one  can  conclude 


2  m+1 
p=l r=l 

where  we  define  Cf  and  f{  j,p)  as  follows  ^  i  j,  m,  p  e  N. 


m>  j  >0 
m<  j 

j  =  0 


f{j,p)  ■=  {ai  +  {-l)-’a2  j>0,p=l 

-bi  +  i-l)tb2  j>0,  p  =  2. 


(11) 


(12) 

(13) 


One  can  see  from  the  right-hand- side  of  (11),  that  the  time  derivative  of  Pm  is  a 
linear  combination  of  the  moments  Pr,  up  to  order  r  =  m  -f  1.  Hence,  if  one  stacks 
all  moments  in  an  infinite  vector  /ico  =  [/ii,  /i2,  ■  ■  ■  ]^>  its  dynamics  can  be  written  as 


/too  =AooAtoc,  (14) 

for  some  infinite  matrix  Aoo.  Let  p  =  [/ii ,  /i2, . . . ,  /i„] ^  G  M"  contains  the  top  n  ele¬ 
ments  of  Poo.  Then,  using  (1 1)  the  evolution  of  p  is  given  by 


p  =Ap  +  Bpn+u 


(15) 


^  n!  denotes  the  faetorial  of  n. 
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for  some  nxn  and  nx  \  matrices  A  and  B  which  have  the  following  structure 


where  the  *  denotes  non-zero  entries.  Our  goal  is  to  approximate  (15)  by  a  finite¬ 
dimensional  nonlinear  ODE  of  the  form 

V  =Av  +  B%+i{v),  v  =  [vi,V2,...,v„]^  (17) 

where  the  map  (p„+i  :  M"  — M  should  be  chosen  so  as  to  keep  v(t)  close  to  /i(t). 
This  procedure  is  commonly  referred  to  as  moment  closure.  We  call  (17)  the  trun¬ 
cated  moment  dynamics  and  (/i)  the  moment  closure  function  for  jln+i- 

When  a  sufficiently  large  but  finite  number  of  derivatives  of  /i(t)  and  v{t)  match 
point-wise,  then,  the  difference  between  solutions  to  (15)  and  (17)  remains  close 
on  a  given  compact  time  interval.  This  follows  from  a  Taylor  series  approximation 
argument.  To  be  more  precise,  for  each  5  >  0  and  integer  N,  there  exists  T  G  M,  for 
which  the  following  result  holds:  Assume  that  for  every  to  >  0, 

At(to)  =  v(to)  and  Vf  e  {1,...,N}  (18) 

where  and  represent  the  time  derivative  of  /i(t)  and  v(t)  along  the 
trajectories  of  system  (14)  and  (17)  respectively.  Then, 

||AtW-v(t)||<5,  VtG[to,r],  (19) 

along  solutions  of  (14)  and  (17),  where  /i  denote  the  first  n  elements  of  In  the 
next  section  we  use  (18)  to  construct  moment  closure  functions  (/i). 


4  Separable  Derivative-Matching  Moment  Closures 

In  this  section  we  construct  truncated  moment  dynamics  (17)  for  the  stochastic  lo¬ 
gistic  model  using  (18).  After  replacing  (15)  and  (17)  in  (18),  equality  (18)  becomes 
a  PDE  on  We  will  seek  for  solutions  cpn+i  to  this  PDE  that  have  the  following 
separable  form 

<+i(m)  =  n  (20) 

m=l 
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for  appropriately  chosen  constants  G  M.  As  we  will  see  in  Section  6.2,  this  sep¬ 
arable  structure  ensures  that  the  steady-state  solutions  for  the  truncated  moment 
dynamics  (17)  are  unique,  positive  real,  and  hence,  biologically  meaningful.  In  the 
sequel  we  refer  to  such  (p^^j(/i)  as  a  separable  derivative-matching  moment  clo¬ 
sure  function  for  jln+l- 


One  can  see  that  the  infinite  vector  G  can  be  expressed  as 


’e[x(0]' 

X 

E[x(t)2 

oo 

=  I 

x^ 

E[x(t)3] 

X=1 

ftW- 


Hence,  the  infinite  vectors  =  [x,x^  ,  which  corresponds  to  a  determin¬ 

istic  distribution,  i.e.  x(t)  =  x  with  probability  one,  form  a  natural  basis  for  Ooo.  In 
particular,  we  will  find  constants  Jm  which  satisfy  (18)  for  each  vector  /ioo(to)  be¬ 
longing  to  this  basis,  i.e.  for  the  class  of  deterministic  initial  conditions.  However, 
often  it  is  not  possible  to  find  Jm  for  which  (18)  holds  exactly.  We  will  therefore 
relax  this  condition  and  simply  demand  the  following 


M(<o)  =  v(,„)  and  =  +  (21) 

Vf  G  {1,2},  where  each  element  of  the  vector  e,(x(to))  is  a  polynomial  in  x(to). 

One  can  think  of  (21)  as  an  approximation  to  (18)  that  is  valid  as  long  as  \t=tQ 
dominates  E[e,-(x(to))]. 


The  following  theorem  summarizes  the  main  result,  the  proof  of  which  is  given  in 
Appendix  B. 


Theorem  2:  Let  7m,  m  G  { 1 , . . . ,  n}  be  chosen  as 


Ym 


(-1) 


n—m 


c 


n-\-l 

m 


(22) 


Then,  for  every  deterministic  initial  condition  v(to)  =  pito)  =  [-^o,  •  •  •  which 
corresponds  to  x(to)  =  xq  with  probability  one,  we  have 


dp{t),  _dv{t), 

dt  dt 

d^p(t)  I  d^v(t)  I  ,  . 


(23a) 

(23b) 


where  and  ^^2^  represent  the  time  derivative  of  p{t)  and  v(t)  along  the 
trajectories  of  the  systems  (14)  and  (17),  respectively,  and  the  n^^  element  of  the 
vector  e2{xo)  is  a  polynomial  in  xq  of  order  2  with  all  other  elements  being  zero.  ■ 
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Remark  1.  Using  (11)  it  can  be  shown  that  is  a  linear  combination  of  mo¬ 

ments  of  X  up  to  order  n  +  2.  Thus,  ^  \t=tQ  is  a  polynomial  in  xq  of  order  n  -f-  2, 
and  hence,  for^  a:o  >>  1,  e^ixo) =  0{xq'^),  where  e^ixo)  is  the 
element  of  e2(-^o)-  Hence,  the  term  dominates  e|(a:o)  by  Xq'^.  ■ 

Remark  2.  It  can  be  verified  that  the  separable  derivative-matching  moment  clo¬ 
sures  also  matches  derivatives  of  order  higher  than  2  in  (23)  with  small  errors.  For 
example  for  n  G  {2, 3, 4}  and  /  G  (2, . . . , 9},  we  have 


d^lit) 

1  dv{t)^ 

dt 

dt 

d‘il{t) 

1  _  ^'v(t)  1 

dd 

II 

1 

o 

II 

+  £i{xo), 


(24a) 

(24b) 


where  the  element  of  e;(vo)  is  a  polynomial  in  xq  of  order  m  —  n  +  i,  for  m  — 
n  -|-  /  >  2  and  equal  to  zero  otherwise.  We  conjecture  that  the  above  equality  holds 
Wn  eN  and  Vz  G  N  but  we  only  verified  it  for  n  up  to  4  and  i  up  to  9.  As  above,  the 

il(t) 

elements  of  vector  '  \t=tQ  dominate  the  corresponding  elements  of  ei{xo)  by  Vq  ”, 
and  hence,  with  increasing  n,  the  truncated  moment  dynamics  v{t)  should  provide 
a  more  accurate  approximations  to  the  lower  order  moments  /i(t).  ■ 


5  Distribution  based  Moment  Closures 


Most  moment  closure  techniques  that  appeared  in  the  literature  start  by  assuming 
a  specific  class  of  distributions  D  for  the  population,  and  use  this  assumption  to 
express  higher  order  moments  as  a  function  of  the  lower  order  ones.  We  refer  to 
such  a  moment  closure  function  as  being  consistent  with  the  distribution  D  and 
define  them  as  follows. 

Definition :  Let  D  be  a  class  of  distributions  parameterized  by  m  parameters  (^i , . . . ,  qm)  G 
with  the  order  moment  given  in  terms  of  the  , . . . ,  as  follows 

=  VfcG  {1,2,...}.  (25) 

The  moment  closure  function  (pj{}^(/i)  for  Hn+i  is  said  to  be  consistent  with  the 
distribution  D  if,  for  every  (^i , . . . ,  q^)  G  one  has  that 


l^n+l  —  fn+lidlj  ■  ■  ■  idm)  — 


(26) 


^  0(.)  denotes  order  of. 
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Table  2 

Unique  Moment  Closure  Funetions  for  the  m*^  order  truneation  [n  =  m)  and  different  dis¬ 
tributions  D. 


D 

m 

Unique  moment  elosure  funetion 

Normal 

2 

(pKix)  =3  H2Ti-2iaf 

Lognormal 

2 

^3iT)  =  ^ 

Pi 

Poisson 

1 

<P2(P)  =P?  +  Pl 

Binomial 

2 

<p|(p)=2^^  {1x2  pf)  +  3AiiP2 

Pi 

where 


Ml 

fiiqu- 

■  ■  ,<?m) 

filial  1  ■ 

■  ■  7  ^m) 

(27) 


For  well  known  elasses  of  distributions  —  sueh  as  lognormal,  normal,  poisson,  or 
binomial  —  we  simply  say  that  is  the  lognormal,  normal,  poisson,  or  binomial 
moment  elosure  funetion. 


5. 1  Techniques  for  obtaining  (p^_^  j 


Generieally,  when  the  dimension  m  of  the  parameter  spaee  ^  is  the  same  as  the 
dimension  n  of  the  domain  of  the  moment  elosure  funetion  (/i),  the  funetional 
equation  (26)-(27)  in  the  “unknown”  (■)  has  a  unique  solution  ^  .  In  faet,  to  de¬ 
termine  (/i)  one  ean  start  by  solving  (27)  for  <71 , . . . ,  in  terms  of  /ii , . . . ,  /i^, 
and  then  substituting  these  baek  in  (26)  to  obtain  a  unique  moment  elosure  fune¬ 
tion  As  we  ehoose  D  to  be  normal  (Whittle,  1957),  log-normal  (Keeling, 

2000),  poisson,  or  binomial  (Nasell,  2003a),  this  proeedure  results  in  the  different 
moment  elosure  funetions  shown  in  Table  2. 

Diffieulties  arise  when  the  dimension  m  of  the  parameter  spaee  ^  is  strietly  smaller 
than  the  dimension  n  of  the  domain  of  the  moment  elosure  funetion  (p^j(/i),  be- 
eause  in  this  ease  the  funetional  equation  (26)-(27)  does  not  have  a  unique  solu¬ 
tion  ^  and  one  ean  find  infinitely  many  moment  elosure  funetions  eonsistent  with 

^  The  reader  is  invited  to  eonvinee  herself  of  this  by  imagining  that  all  funetions  are  loeally 
linear,  in  whieh  ease  (p^j(/r)  is  represented  by  a  veetor  with  n  unknowns.  Sinee  (26)  must 
hold  for  a  (loeal)  basis  of  we  have  exaetly  m  equations  to  determine  the  n  =  m  unknowns. 
^  In  eontrasts  to  the  previous  ease,  we  now  have  n  unknowns  to  represent  the  loeal  lin- 
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the  same  family  of  distributions.  However,  there  is  a  strong  incentive  to  consider 
this  case  because,  as  well  shall  see  shortly,  large  values  for  n  generally  lead  to  sig¬ 
nificantly  more  accurate  moment  closures.  In  the  sequel,  we  illustrate  some  options 
for  moment  closures  with  n>  m. 

Example  1:  Consider  the  class  of  poisson  distributions  characterized  by  their  ex¬ 
pected  value  6  (m=  1).  Their  moments  are  given  by 

Ati=/i(0):=0 

At2  =  /2(0):=  0(1  +  0) 

At3=/3(0):=  0(1+30  +  02),... 

Nasell  (2003a)  proposed  the  following  poisson  moment  closure  function  for  n  =  2: 

(pP\ll)=lli  +  3piP2-2lll,  (28) 


for  which  it  is  straightforward  to  verify  that  (26)-(27)  holds  because 

P3  =  e{i  +  3e  +  e^)  =  (pP\p),  p  =  [0 ,  e{i  +  e)f. 

However,  an  alternative  choice  for  the  poisson  moment  closure  function  that  also 
satisfies  (26)-(27)  is  given  by 


(pfilii)  =ji2-  ^2  +  3jiiji2  -  2Atf ,  (29) 

which,  as  we  will  see  in  the  next  section,  performs  better  that  (28).  The  explanation 
for  this  lies  in  the  fact  that  (29)  has  better  derivative  matching  properties  than  (28), 
in  the  sense  of  (21).  In  the  sequel  we  refer  to  (28)  and  (29)  as  the  Nasell-poisson 
and  new-poisson  moment  closure  function,  respectively. 

Example  2:  Consider  now  the  class  of  normal  distributions  parameterized  by  their 
mean  co  and  variance  (m  =  2).  Their  moments  are  given  by 


Ail  =/i(0,c7) 
Ii2  =  f2{0,a) 

P4  =  f4{9,a) 


=  CO 

=  ©2  -P  cy2 

=  co{co^  +  3a^) 

=  +  6co^a^  +  3a‘^, . . . 


For  n  =  3,  any  function  of  the  following  form  is  a  normal  moment  closure  function 
for  /i4 


(p4{p)-=  4AtiAt3  +  3pl  -  np2lil  +  6Atf  Eh{pi,p2-  , Ms  -  3p2lii  + 


earization  of  (/r)  but  a  (local)  basis  of  ^  only  provides  m  <n  equations. 
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where  h{x,y,z)  is  any  function  with  the  property  that  h{x,y,0)  =  0.  To  verify  that 
this  is  so,  we  note  that  (26)-(27)  holds  because  for 


/i  =  [to,  (0^  +  0^,  (o{(o^ +  3a'^)Y : 

we  obtain  /z(/ii,/i2  — —  3/i2/ii  +2/i^)  =  h{co,  C7^,0)  =  0  and  therefore 

(P4(/i)  =  +  6C0^G^  +  3(7  =  /i4. 


Example  3:  Finally  consider  the  class  of  lognormal  distributions  characterized  by 
the  parameters  a  >  0,  j8  >  0  (m  =  2),  whose  moments  are  given  by 

Pk  =  :=  e  {1,2, ■  ■  ■ }.  (30) 

For  every  n  >2,  the  separable  derivative-matching  moment  closure  functions  de¬ 
fined  by  (20),  with  coefficient  given  by  (22)  in  Theorem  2  are  lognormal  moment 
closure  functions  for  Pn+i-  We  can  verify  this  by  noting  that  (26)-(27)  holds  be¬ 
cause,  from  (30)  and  (B.l)  used  in  the  proof  of  Theorem  2  (Appendix  B)  we  con¬ 
clude  that 


=  aiK=irn,cT)p{r„=irm{2C^+cT})^ 

where  we  used  the  facts  that  k  =  C\,k^  =  2C\  +  C\,\/k  eN. 


5.2  Cumulant  closure  functions 


Most  of  the  literature  on  moment  closure  prefers  to  work  with  a  vector  K  =[ki,  . . .  ,Kn\ 
where  Kn{t)  is  the  n^^  order  cumulant^  ,  instead  of  the  previously  introduced  vector 
p  of  uncentered  moments  in  (15).  Then,  instead  of  doing  moment  closure  one  per¬ 
forms  cumulant  closure  by  approximating  Kn^i  by  a  nonlinear  function 
of  Ki, . . . ,  tCn,  which  we  refer  to  as  the  cumulant  closure  function.  The  disadvan¬ 
tage  of  working  with  K  instead  of  p  is  that  the  dynamics  of  K  is  always  nonlinear. 
However,  for  ease  of  comparison  with  other  papers,  we  provide  in  Table  2  the  cu¬ 
mulant  closure  functions  corresponding  to  the  different  moment  closure  functions 

2  The  order  eumulant,  Kn  is  given  as  follows  in  terms  of  the  uneentered  moments 

Ki=pi,  K2=p2-pf 

K2=  P3-3piP2+2pf,  K4=p4-4piP3-3pi  +  l2p2Pi-6p^,... 
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Table  3 

Moment  Closure  Funetions  (MCF)  for  seeond  order  truneation  {n  =  2)  and  eorresponding 
Cumulant  Closure  Funetions  (CCF)  for  /ta  and  K^,  respeetively,  eorresponding  to  the  differ¬ 
ent  Moment  Closure  Teehniques  (MCT)  diseussed  in  this  paper.  SDM  refers  to  separable 
derivative-matehing . 


MCT 

SDM 

Normal 

Lognormal 

Nasell-Poisson 

New-Poisson 

Binomial 


MCF 

•PsCm)  =  ^ 

Ml 

(pI{Ix)  =  31X21^1 -2ixf 


(pP\lx)  =  IXi+3iXiiX2-2ixf 

(pf{lx)  =  H2-l^i+3iXiiX2-2ix^i 

-  (Ai2  -  At?)  +  2 Hi  1X2  -  2ixf 
_ Ml _ 


CCF 


2  T~ 


■CK 


^I{k)=0 


(j)P\K)  =  Ki 

^f(K)  =  K2 

'/’I (M)  =  2^  -  K-2 

_ M _ 


discussed  so  far  for  n  =  2.  We  use  superscripts  s,  I,  g,  pi,  p2  and  b  to  denote  sep¬ 
arable  derivative-matching,  lognormal,  normal,  Nasell-poisson,  new-poisson  and 
binomial  moment  closure  functions,  respectively. 


6  Comparison  of  Moment  closures 


In  this  section,  we  introduce  two  criteria  to  compare  the  different  moment  closure 
techniques.  The  first  criterion  is  the  error 


•  Mn+l(0  ^)i+i(M(0) 


oo 


I 


— ^e«+i 


where 


■  rfAWiWi  ^'<Pn+i(M(0)| 

n+lK-^O)  ■ 


(31) 


(32) 


We  call  e^_|_i(.t:o)  the  derivative  matching  error.  Ideally,  one  would  like  to  have 
^M+i  (-^o)  =  but  as  already  pointed  out  in  Section  4,  this  is  generally  not  possible. 
With  deterministic  initial  conditions  /foo(to)  =  Theorem  2,  the 

derivative  matching  error  is  typically  a  polynomial  in  xq.  The  lesser  the  order  of 
this  polynomial,  the  better  is  (p„_|_i(/i)  in  approximating  /i„+i. 


The  second  criterion  is  the  steady-state  solution  of  the  truncated  moment  dynamics 
(17).  In  particular,  we  are  interested  in  determining  if  there  exists  a  unique  non¬ 
trivial  positive  real  steady-state  which  is  physically  meaningfull.  This  is  important. 


14 


because  it  is  well-known  that  normal  moment  closures  can  have  spurious,  imagi¬ 
nary  and  sometimes  even  stable  negative  steady-states,  which  lead  to  biologically 
meaningless  solutions  (Keeling,  2000). 

6. 1  Derivative  Matching  Error 

6.1.1  Moment  Closures  for  n  =  2 

We  recall  from  Table  3  that  and  therefore  we  do  not  need  to  dis¬ 

cuss  lognormal  moment  closure  separately.  By  substituting  cp^in), 

(Pg  ^(/i)  and  from  Table  3  in  (31)-(32),  one  obtains  the  corresponding  deriva¬ 
tive  matching  errors,  which  will  be  denoted  using  the  appropriate  superscripts. 

Using  Table  3  and  symbolic  manipulation  in  Mathematica,  we  can  show  that 


^eO(.ro)  =  ^e^ixo)  =  Ph^{xo)  =  ^s^(xo)  =  0  (33a) 

^^£30(^0)  =  -^0-  (33b) 

*£3(.ro)  e^xo(i  +  l),  *  ^  {^,g,pi,p2,h},  Vie  {1,2,...}  (33c) 

where  Pxo(J)  denotes  the  set  of  polynomials  in  xq  of  order  j.  Since  ^^£3  (.ro)  =  —xq. 


the  Nasell-poisson  moment  closure  function  will  have  a  large  initial  error,  espe¬ 
cially  for  large  initial  conditions,  when  compared  to  all  other  moment  closure  func¬ 
tions.  For  all,  /  G  (1,2, . . .},  all  of  these  moment  closure  functions  match  deriva¬ 
tives,  with  the  derivative  matching  error  being  of  the  same  order  in  xq.  The  simu¬ 
lation  results  discussed  below  show  that  with  the  exception  of  Nasell-poisson  mo¬ 
ment  closure  function,  which  consistently  provides  the  worst  estimates,  all  other 
moment  closure  functions  perform  fairly  well. 

Example:  We  consider  the  stochastic  logistic  model  with 

ai  =  .30,  a2  =  .02,  bi  =  .015,  62  =  .001,  (34) 

which  is  used  by  Matis  et  al.  (1998)  to  model  the  population  dynamics  of  the 
African  Honey  Bee.  Using  (17)  with  the  matrices  A  and  B  computed  in  (11),  we 
have  the  following  truncated  moment  dynamics 

Vi  0.28  -0.016  Vi  0 

+  <P3(V). 

V2  .32  .546  V2  -0.032 

The  time  evolution  of  the  moments  corresponding  to  different  moment  closure  tech¬ 
niques  is  obtained  by  substituting  the  appropriate  moment  closure  function  from 
Table  3  in  place  of  <P3(v).  Figure  1  plots  the  different  errors  (31)  during  the  time 
interval  [0,2.5]  for  jcq  =  5.  The  error  is  approximated  by  the  first  nine  terms  of 
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Fig.  1.  Evolutions  of  the  error  e„+i  in  (31)  for  the  different  moment  elosure  teehniques  in 
Table  3,  with  parameters  as  in  (34)  and  xq  =  5. 


the  series  in  (31).  The  letters  s,  g,  pi,  p2  and  b  are  used  to  denote  the  errors  eor- 
responding  to  separable  derivative-matehing,  normal,  Nasell-poisson,  new-poisson 
and  binomial  moment  elosure  functions,  respectively.  In  order  to  evaluate  the  per¬ 
formance  of  these  moment  closures  techniques  past  the  initial  period,  we  also  com¬ 
pute  the  exact  evolution  of  the  moments.  This  is  only  possible  because  the  popu¬ 
lation  limit  t/  =  25  is  small  and  one  can  obtain  the  exact  solution  by  numerically 
solving  the  Kolmogorov  equation  (6).  Figure  2  contains  plots  of  the  mean  and  vari¬ 
ance  errors,  respectively,  for  the  different  moment  closure  functions  with  v:o  =  5  and 
xq  =  20.  The  letters  s,  g,  pi,  p2  and  b  are  used  to  denote  the  errors  corresponding  to 
separable  derivative-matching,  normal,  Nasell-poisson,  new-poisson  and  binomial 
moment  closure  functions,  respectively.  For  xq  =  20  the  binomial  moment  closure 
function  provides  the  best  estimate  both  initially  and  at  steady-state,  whereas  for 
v:o  =  5  the  new-poisson  moment  closure  function  does  best  initially,  but  the  bino¬ 
mial  moment  closure  function  continues  to  provide  the  most  accurate  steady-state 
estimate.  As  one  would  expect  from  (33),  the  Nasell-poisson  moment  closure  func¬ 
tion  performs  the  worst. 


6.1.2  Moment  Closures  for  n>3 

All  distributions  D  discussed  in  Section  5  are  characterized  by  at  most  2  parame¬ 
ters  (m  <  2),  hence  the  corresponding  distribution-based  moment  closure  functions 
are  not  unique  for  n  >  3.  In  this  section  we  illustrate  for  lognormal  and  normal 
moment  closures,  how  derivative  matching  can  be  used  as  a  tool  for  gauging  their 
performances  and  choosing  the  ones  that  yield  the  least  derivative  matching  errors. 
We  also  propose  a  new  moment  closure  function  for  n  =  3,  which  as  we  will  see 
guarantees  better  approximation  at  least  locally  in  time,  as  compared  to  other  ones. 
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(c)  variance  error  for  xq  =  5  (d)  variance  error  for  xq  =  20 


Fig.  2.  Evolutions  of  the  mean  error  /Ti  —  Vi  and  of  the  variance  error  (/r2  —  fij^)  —  (V2  —  Vj^) 
for  the  different  moment  closures  in  Table  3  for  n  =  2,  with  parameters  as  in  (34)  and  xq  =  5 
(left)  xo  =  20  (right). 

Lognormal  moment  closure:  One  can  see  from  Example  3  of  Section  5,  that  a 
family  of  lognormal  moment  closure  functions  for  is  given  by 

=  £c;'(„  =  cf‘.  £c%  =  cf‘.  (35) 

m=l  m=l 

The  derivative  matching  errors  for  these  moment  closure  functions  are  given  by 
^^n+li^o)=^i  ^^n+li^o)  ^  Px+oifl  +  i+^~w)j  V/g{1,2,...}  (36) 

where  the  constant  w  G  { 1 , 2, . . . ,  n}  is  such  that 

m=\ 

VfcG{w+l,...,n} 

m=i 

with  higher  values  of  w  correspond  to  better  approximations.  From  (35)  and  (36), 
we  conclude  that  the  lognormal  moment  closure  guarantees  that  w  >  2.  But  from 
(B.l)  we  see  that  this  can  be  improved  and  we  actually  set  w  =  n  for  the  separable 
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derivatives -matching  moment  closure,  for  which 

^e„V(^o)  =  0,  ^ei+,ixo)eP,+o{i  +  l).  Vi  e  {1,2,...}.  (37) 

Thus,  the  separable  derivative-matching  moment  closure  function  leads  to  the  best 
estimate  for  jJ-n+i,  in  the  sense  that  it  has  the  least  derivative  matching  error  among 
all  the  moment  closure  functions  consistent  with  the  lognormal  distribution. 


Normal  moment  closure:  We  first  restrict  our  attention  to  the  case  n  =  3.  We 
recall  from  Example  2  that  a  family  of  normal  moment  closure  functions  for  114 
were  given  by 

4^1  At3  +  3At2  -  12At2Ati^  +  +  m  -  +  2Ati ) • 

Its  straight  forward  to  see  that  the  function  h(K'i,  K^,  K‘3)  corresponds  to  the  cu- 
mulant  closer  function  for  K4.  However,  one  needs  to  be  careful  in  picking  the 
functions  h.  To  demonstrate  this  we  consider  the  following  normal  moment  closure 
functions 


<P|(At)  =  4AtiAt3  +  3At2  -  12;U2At?  +  (38) 

(m)  =  4Ati  At3  +  3At|  -  12At2At?  +  6/1^  +  P3-  (39) 

(pfiP)  =  4AtiAt3  +  3At|  -  12At2At?  +  -h  {113  -  3AtiAt2  +  2Atf)Ati  (40) 

which  correspond  to  cumulant  closure  functions  h(K‘i,  K^,  K3)  =  0,  K3  and  K‘3K‘i 

respectively.  Using  symbolic  manipulation  in  Mathematica  we  have  the  following 
derivative  matching  errors  for  the  normal  moment  closure  functions  (38)-(40): 

*e}(xo)  e /!co(f-l- 1),  *  =  {g,gl} 
s^eiixo)eP,,ii  +  2) 

for  all  /  G  (1, 2, . . .}  and  zero  for  i  =  0.  Thus,  the  normal  moment  closure  functions 
(38)-(39)  have  derivative  matching  errors  of  the  same  form  as  (37),  and  provide 
better  approximates  for  as  compared  to  (40),  for  which  the  derivative  matching 
errors  are  an  order  higher. 

In  general  for  n  >  3,  a  family  of  normal  moment  closure  functions  for  is 
obtained  using  the  cumulant  closure  function 

=h(Ki,K-2,...,K-„) 

where  the  function  h  is  zero  if  any  of  K3,. . .  ,Kn  is  zero.  Using  symbolic  manipula¬ 
tion  in  Mathematica  we  found  that  choosing  cumulant  closure  functions  as 

n 

h(Ki,K-2,...,K-„)  = 

d=3 
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for  some  constants  /^,  typically  lead  to  derivative  matching  error  of  the  form  (37) 
and  provides  good  estimates  for 


Zero  first-order  error  moment  closure:  We  now  propose  a  slight  modification 
of  the  normal  moment  closure  function  (38)  which  is  given  by 

(plin)  =  AiiiHs  +  3ju|  -  12jU2Ati^  +  +  M2  -  At?  (41) 

and  corresponds  to  the  cumulant  closure  function  K2  for  K4.  It  yields  the  following 
derivative  matching  errors 

^£4(^:0)  =0,  /g{0,  1}  ^£\{xq)  &  Px^ii+l),  Vi  >  2. 

Hence,  unlike  the  separable  derivative-matching  and  normal  moment  closure  func¬ 
tions,  (plin)  yields  zero  and  order  derivative  matching  errors,  and  hence, 
provides  the  best  estimates  for  as  compared  to  them,  at  least  near  t  =  0. 


In  order  to  gauge  the  performances  of  the  above  moment  closure  functions  we 
consider  the  stochastic  logistic  model  with  parameters  as  in  (34),  n  =  3  and  xq  =  20. 
We  recall  from  Table  1  that  for  n  =  3,  we  have  the  following  separable  derivative¬ 
matching  moment  closure  function 

(p'liP)  =  (42) 

M2 

Using  (1 1),  we  have  the  following  truncated  moment  dynamics 


Vi 

0.28 

-0.016 

0 

Vi 

0 

V2 

= 

.32 

.546 

-0.032 

V2 

+ 

0 

V3 

.28 

.944 

.798 

V3 

-0.048 

<P3(V). 


Substituting  the  moment  closer  functions  (38)-(42)  in  place  of  <P3(v),  we  obtain 
the  corresponding  approximate  time  evolution  of  moments.  Figure  3  contains  plots 
of  the  mean,  variance,  and  third  cumulant  errors.  The  letters  g,  gl,  g2,  z  and  s 
represent  these  errors  for  the  moment  closure  functions  (38)-(42),  respectively.  As 
expected  the  normal  moment  closure  functions  (38)-(39)  perform  much  better  than 
the  normal  moment  closure  function  (40).  The  separable  derivative-matching  mo¬ 
ment  closure  function  (42)  also  provide  good  estimates.  One  can  also  see  that  the 
zero  first-order  error  moment  closure  function  (41),  which  guarantees  the  best  ap¬ 
proximation  near  t  =  0  actually  provides  in  this  case  the  most  accurate  estimate  for 
/i4  for  all  time.  As  can  be  seen  from  Figure  2  and  3,  the  mean  and  variance  errors 
for  Xq  =  20  with  n  =  3  are  an  order  of  magnitude  smaller  as  compared  to  the  ones 
with  n  =  2. 
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(a)  mean  error  for  xq  =  20 


t 

(b)  varianee  error  for  xq  =  20 


(e)  third  eumulant  error  for  xq  =  20 


Fig.  3.  Evolutions  of  the  mean  error  /Ti  —  Vi,  varianee  error  (1^2  — l-lf)  —  (V2  — vf),  and  third 
eumulant  error  (/ra  —  3/X2fii  +  2/rj^)  —  (V3  —  3V2Vi  +  2Vj^)  for  the  different  moment  elosure 
funetions  (38)-(42)  for  n  =  3,  with  parameters  as  in  (34)  and  xq  =  20. 

6.2  Steady-State  Solutions  of  the  Truncated  Moment  Dynamics 


We  now  look  at  the  steady-state  solutions  of  the  truncated  moment  dynamics  (17) 
for  the  different  moment  closure  techniques. 


6. 2. 1  Separable  Derivative-Matching  Moment  closure 

Consider  the  truncated  moment  dynamics  of  order  n  G  N>2  with  moment  closure 
functions  as  given  in  Table  1.  From  (17),  at  steady-state  we  have 

0=Av^'H+5(p^^j(v"(oo))  (43) 

where  V^(oo)  denotes  the  steady-state  solution.  Using  (16)  at  steady-state  we  have 

v|(oo)  =  ciVi(oo)  (44a) 

(44b) 
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(44c) 


<P«+l(v"(oo))  =c„V[(oo) 


for  some  positive  real  numbers  ci, . . .  ,c„.  Using  the  above  equalities  and  Table  1 
we  eonelude  that 


'p«v'H)  =  c5,  «>«vh)  = 'Pl(v'W)  = 


„10„5 

^10 


(45) 


for  n  G  {2,3,4,...}  respeetively.  Substituting  (45)  in  (44)  yields  the  following 
unique  non-trivial  solution  for  v{(oo) 


C2 


v;(~)  =  ■!  # 


„10_5 


n  =  2 
n  =  3 
n  =  4 


and  the  eorresponding  V2(°o), . . . ,  v^(oo)  ean  be  ealeulated  from  (44).  Henee,  the 
separable  derivative-matehing  moment  elosure  always  yields  a  unique  non-trivial 
positive  real  steady-state  Vn  G  N>2.  In  terms  of  the  parameters  ai,  bi,  a2  and  b2, 
the  eonstants  ci,  C2  and  C3  are  given  as  follows 


Cl  =K,  C2  =  K^  +  0^, 

ai-a2  2_^1^2  +  ^1«2 

^1+^2’  ^  ibi+b2Y 


C2,=K^  +  3Ko^  -1-  00^ 
_  _  b2-bi 
^  b2  +  bi' 


and  henee, 


V{(oo) 


K 


1  + 


<y^ 


n  =  2 

n  =  3. 


6.2.2  Other  Moment  closure  techniques 

In  this  seetion  we  will  see  that  the  moment  elosure  funetions  whieh  laek  the  “sep¬ 
arable  strueture”  of  (20),  ean  lead  to  seenarios  of  biologieally  meaningless  steady- 
state  solutions  for  the  truneated  moment  dynamies  (17).  Substituting  the  moment 
elosure  funetions  from  Table  2  in  (43),  we  get  the  following  non-trivial  steady-state 
solutions  for  n  =  2  : 


vf(oo)  =i^ 
vf^(oo)  =i^ 


-3 

1 

/ 

±-l 

1- 

.4 

4 

-3 

1 

/ 

±-l 

1- 

.4 

4 
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8a^\  2' 
~W)  . 

8(a2-l)x 

) 


1 

2 


pi^-1  ^  1 
-  4  ^4 


Sct^n  i- 

~w) . 


From  the  above  steady-states  we  conclude  the  following  for  n  =  2: 

(1)  The  binomial  moment  closure  function  leads  to  a  unique  non-trivial  attracting 
real  steady-state,  which  can  be  negative  for  a  range  of  parameters. 

(2)  Normal,  Nasell-poison  and  new-poison  moment  closure  functions,  yield  two 
non-trivial  steady-states,  with  the  one  with  the  —  sign  being  a  “spurious  steady- 
state”.  Following  the  definition  by  Nasell  (2003b),  a  steady  state  is  “spuri¬ 
ous”,  if  limM^oo  Vi  (oo)  7^  K,  where  and  K  are  both  0{M).  For  n  =  2,  all 
these  “spurious  steady-states”  happen  to  be  unstable,  and  hence,  the  truncated 
model  will  not  converge  to  them. 

When  the  parameters  are  chosen,  such  that  the  term  under  the  square  root 
sign  is  negative,  then  both  the  non-trivial  steady-states  would  be  imaginary, 
and  hence,  biologically  meaningless. 

In  general  for  n  >  3,  the  normal  moment  closure  functions  introduced  in  Sec¬ 
tion  6.1.2  yield  K>3  non-trivial  steady-states  with  1  of  them  being  “spurious”. 
For  example,  for  n  =  3,  the  normal  and  the  zero  first-order  error  moment  closure 
functions  (38)  and  (41)  yield  3  non-trivial  steady-state  solutions  given  as  the  roots 
of  the  third  order  polynomials 

(3cf  -f  4c2)vi(oo)  -  12ciVi(oo)^-h6vi(oo)^  =  C3, 


and 


{3cI  +  4c2  -  l)vi(oo)  -  12ciVi(oo)^-h6vi(oo)^  =  C3  -ci, 

respectively,  with  2  of  them  being  “spurious”.  Also  for  a  range  of  parameter  val¬ 
ues,  biologically  meaningless  scenarios  of  a  combination  of  imaginary  and  negative 
steady  states  can  happen.  On  the  other  hand,  the  separable  derivative-matching  mo¬ 
ment  closure  function  with  a  unique  non-trivial  positive  real  steady-state  Vn  G  N>2 
has  a  clear  advantage. 


7  Conclusion  and  Future  Work 


A  procedure  for  constructing  moment  closures  for  the  stochastic  logistic  model 
was  presented.  This  was  done  by  first  assuming  a  separable  form  for  the  moment 
closure  function  (p„^i{v),  and  then,  matching  its  time  derivatives  with  /i„+i,  at 
some  initial  time  to-  We  showed  that  for  initial  conditions  x(to)  =  xq  with  prob¬ 
ability  one,  there  exists  a  unique  separable  derivative-matching  moment  closure 
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function  for  which  the  derivative  matching  error  is  a  polynomial  in  xq  of  order 
z  +  1  for  all  z  G  {1,2, . . .}  and  zero  for  z  =  0.  Explicit  formulas  to  construct  these 
moment  closure  functions  were  provided.  Comparisons  with  alternative  moment 
closure  techniques  available  in  literature  were  carried  out  based  on  transient  per¬ 
formance  and  steady-state  solutions  of  the  truncated  moment  dynamics,  which  led 
to  the  following  conclusions: 

(1)  Derivative  matching  can  be  used  as  a  effective  tool  for  gauging  the  perfor¬ 
mance  of  moment  closure  functions.  We  showed  that  for  n  =  2,  with  the  ex¬ 
ception  of  the  Nasell-poisson,  all  other  moment  closure  functions  in  Table  2 
perform  derivative  matching  and  give  fairly  good  estimates  of  /is.  For  n>  3, 
we  showed  that  among  the  moment  closure  functions  consistent  with  the  log¬ 
normal  distribution,  the  separable  derivatives-matching  moment  closure  func¬ 
tion  provides  the  best  estimate  for  /i„+i.  We  also  proposed  families  of  normal 
moment  closure  functions  that  perform  derivative  matching  and  result  in  good 
estimates.  For  zz  =  3,  a  new  zero  first-order  error  moment  closure  function  was 
also  proposed  which  guaranteed  better  approximations,  at  least  locally  in  time, 
as  compared  to  the  other  moment  closure  functions  discussed  in  this  paper. 

(2)  The  separable  derivative-matching  moment  closure,  always  yields  a  unique 
non-trivial  positive  real  steady-state  Vzz  G  N>2,  and  hence,  in  some  sense  su¬ 
perior  to  the  other  moment  closures,  which  can  have  spurious,  imaginary  and 
even  stable  negative  steady-states. 

Possible  directions  for  future  research  include  the  development  of  a  systematic  pro¬ 
cedure  to  construct  moment  closure  functions  of  the  form  (41)  which  can  yield  zero 
derivative  matching  errors  up  to  any  order  and  the  extension  of  the  results  in  this 
paper  to  multi-specie  birth-death  Markov  processes.  Primary  results  regarding  the 
latter  can  be  found  in  Hespanha  and  Singh  (2005). 


A  Appendix 


Substituting  (4),  (7)  and  (8)  in  (10)  and  doing  a  binomial  expansion  we  have 

d  j^fj 


dt 


=  E  [(x  +  1)“  -  X™)!]  (x)  -f  (x  -  1)"*  -  x^)x{x)] 


E 


E 


£  CJ[{aix-bix^)  +  {-\y{a2X  +  b2X^)]x"^-^ 
.7=1 

2  m 

£  X;C7/(7,p)x-^+i^ 

_p=lj=\ 


(A.l) 


where  CJ  and  f{j,p)  are  defined  in  (12)  and  (13).  Using  (A.l),  (12)  and  (13)  we 
conclude  that  the  evolution  of  /!„,  Vzzz  G  N,  can  be  written  as 
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2  m 

Am=  ^ 

P=ij=i 
2  m+p~l 

^m+p—  J{m  +  p-r,p)jlr,  r^m-j  +  p 

p=i  r=p 
2  m+1 

=  L  E  CZ+p^rf{m  +  p-pp)Pr-  (A.2) 

p=l  r=l 


B  Appendix 


Proof  of  Theorem  2: 

STEP  1 :  We  first  show  that  with  7^  chosen  as  the  solution  to  the  linear  system  of 
equations 


Cr'=ErmCr,  (B.l) 

m=\ 

equalities  (23)  hold.  Using  (15),  (16)  and  (17)  one  can  see  that  the  following  equal¬ 
ities  imply  (23): 


At«+l(?o)  =  <Pn+l(v(to)) 

4U+lW|  _  ^<Pn+l(v(0)|  /  X 

dt  dt 

where  (.ro)  is  a  polynomial  in  xq  of  order  2. 


(B.2) 

(B.3) 


Let  Ym  be  chosen  as  the  solution  of  (B.l).  We  show  next  that  equalities  (B.2)  and 
(B.3)  hold.  With  initial  conditions  starting  on  the  set  of  deterministic  distributions, 
we  have  /im(to)  =  -^0  •  From  /i(to)  =  v(to),  (20)  and  (B.l)  we  have 


which  proves  (B.2).  From  (A.2)  we  have 

2  n+2 

A«+lW=  E  E^«+l+p-r/('^+^+^“0P)Pr 

p=l r=l 
2  n+2 

fln+M  =  E  E^«+i+p-r/('^+^+^“OPK  +  ie„Vi(^o),  (B.4) 

p=l r=3 

where  (.ro)  is  a  polynomial  in  xq  of  order  2.  Using  (A.2)  and  (17)  we  have 
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2  m+1 


^m+p~  J{m  +  p-r,p)Vr,  V„+i  =  (p^+i(v). 

p=l  r=\ 


(B.5) 


Thus,  from  (20),  (B.l),  (B.5)  and  v(?o)  =  [.^o,  ■  ■  •  5-^0] 


nYT 


^‘Pn+l(v(0) 


dt 


'm 

m=l 
2  n  m+1 

=  £  £  £  Tm vf  ■  ■  ■  1  ■  ■  ■  vJ"C™+p_,/(m  +  ;?  -  r, ;?)  V, 

p=l m=l r=l 

I  =V  v"v  ,  f  +  n  -  r 

lf=fo 


dt 


=  EEE  Tra^m+p-  J{m  +  p-pp)xQ 

p=\ m=l r=l 
2  n  m+1 


=  EEE  7mC„_l_p_  rfim  +  p-pp)^Q 

p=\ m=l r=l 


.n+l+r—m 


With  a  change  of  variable  ^  =  n  +  l  +  r  —  m,  the  last  equality  becomes 


dtn+liyi.t))  ^ 
dt 


2  fj  m+2 

EE  E  rmC+i+p„/(n  +  l+;?-^,;?)+^ 

p=  1  m=  1  q=n+2—m 


2  n+2  n 

^  E  E  E  y^^n+i+p-  gf{n  +  l+p-q,  p)xl  +  |e„Vi  (-^o) , 

p=l  q=3m=l 


(B.6) 


where  2£,!+i(-^o)  is  a  polynomial  in  +0  of  order  2.  Using  (B.l),  equality  (B.6)  re¬ 
duces  to 


^<P«+i(v(0) 

dt 


2  n+2 

E  £  ^n+l+p-^/(^  +  1  +p 


p=l^=3 


^,.P)^0  +  2enVl(-^0)- 


(B.7) 


Comparing  (B.7)  with  (B.4)  one  can  see  that  (B.3)  holds. 


STEP  2:  We  now  show  that  solution  to  (B.l)  is  unique  and  given  by  (22).  Towards 
that  end,  we  observe  that  for  all  z  G  M,  one  can  write  using  binomial  expansion. 


n+l 

ii-(i+z)r‘  =  Ecr‘(-ir(i++ 

V- 0 

n+l  V 

=  l-hEcr^(-l)''ECz"  (B.8) 

v=l  w=0 

Equating  coefficients  for  z'^,  5  G  {1, ...,«}  on  both  sides  of  (B.8)  we  have 


n+l 

o=Ec+'(-i)''cr 


V=1 
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(B.9) 


V=1 

Comparing  (B.9)  with  (B.l)  one  can  see  that  a  solution  to  (B.l)  will  be  (22).  Also 
the  system  of  n  linear  equations  (B.l)  can  be  put  into  the  form 


c  =  n7 


where  7=  [71, . . . , 7„]^,  C  =  and 


n  = 


/^1  /^2 
*^1  •" 

0  cl ...  Cl 


0  0  ...C" 


As  the  upper  triangular  matrix  fl  is  non-singular,  the  above  solution  to  fm  is  unique. ■ 
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